#R version 4.2.1 (2022-06-23)
#Platform: x86_64-w64-mingw32/x64 (64-bit)
#Running under: Windows 10 x64 (build 19044)

library(dplyr) #version 1.0.9
library(readr) #version 2.1.2
library(stringr) #version 1.4.0
library(zoo) #version 1.8-10
library(sandwich) #version 3.0-2
library(lmtest) #version 0.9-40
library(prediction) #version 0.3.14
library(ggpubr) #version 0.4.0
library(estimatr) #version 1.0.0
library(margins) #version 0.3.26
library(broom) #version 1.0.0
library(scales) #version 1.2.0

`%notin%` <- function(x,y) !(x %in% y) 

#GWF predicted plots Figure 1#####
#institutionalization#
d1 <- read_csv("data/gwf_fa_pi.csv") 
d2 <- read_csv("data/gwf_irt.csv") %>% 
  rename(gwf_casename = case) %>%
  select(gwf_casename, transition, gwf_irt) 

gwf <- read_csv("data/panel.csv") %>% 
  left_join(d1, by = c("gwf_casename", "year")) %>% 
  left_join(d2, by = c("gwf_casename")) %>% 
  group_by(gwf_casename) %>% 
  mutate_at(.vars = vars(gwf_incumbent_pi), ~na.locf0(.)) %>% 
  group_by(countryname) %>% 
  mutate_at(.vars = vars(lag_milexp, lgupop, e_migdppcln, e_regionpol, e_total_resources_income_pc, islam50, e_peaveduc, ginibestestimateextrapolatednly, LagRuralInequality, larea, gwf_priorduration, e_coups, EFindex, theta_mean), ~na.approx(., maxgap = 5, rule = 2)) %>% 
  mutate_at(.vars = vars(lag_milexp, lgupop, e_migdppcln, e_regionpol, e_total_resources_income_pc, islam50, e_peaveduc, ginibestestimateextrapolatednly, LagRuralInequality, larea, gwf_priorduration, e_coups, EFindex, theta_mean), ~na.locf0(.)) %>% 
  mutate(v2xps_party = ifelse(is.na(v2xps_party), 0, v2xps_party)) %>% 
  mutate_if(is.numeric, ~replace(., is.nan(.), NA)) %>% 
  mutate(coups = cumsum(e_coups)) %>% 
  select(countryname, year, gwf_casename, gwf_regimetype, v2xps_party, v2pssunpar, gwf_incumbent_pi, gwf_irt, transition,
         e_regionpol, e_migdppcln, lgmilex_cap, lag_milexp, e_peaveduc, lgupop,
         ginibestestimateextrapolatednly, LagRuralInequality, e_total_resources_income_pc, islam50, larea, gwf_priorduration, coups, gwf_number, EFindex, theta_mean) %>%
  mutate(gwf_incumbent_pi = ifelse(is.na(gwf_incumbent_pi), v2xps_party, gwf_incumbent_pi)) %>% 
  rename(incumbent_pi = gwf_incumbent_pi) %>% 
  filter(transition == year, gwf_casename %notin% c("Bolivia 43-46", "Thailand 57-73", "Pakistan 58-71", "Haiti 88-90",
                                                    "Sudan 58-64", "Ecuador 63-66", "Venezuela 48-58")) %>% 
  ungroup() %>% 
  unique()

lm(gwf_irt ~ incumbent_pi + ginibestestimateextrapolatednly + LagRuralInequality + lag_milexp + larea + lgupop + e_migdppcln + e_total_resources_income_pc + islam50 + e_regionpol + e_peaveduc + gwf_priorduration + coups + gwf_number + EFindex + theta_mean, data = gwf) %>% 
  prediction(.) %>% 
  select(incumbent_pi, fitted, se.fitted) %>% 
  ggplot(., aes(x = incumbent_pi, y = fitted)) +
  geom_smooth(color = "black") +
  geom_rug(sides = "b") +
  theme_bw() +
  ylab("Bounded Democracy") +
  ylim(0,1) +
  xlim(0,1) +
  xlab("Party Institutionalization") +
  theme(legend.position="bottom") 

ggsave(filename = "figures/appendix_a_fig1_1.jpg", dpi = 500, width = 5, height = 5) 

#strength models #
d1 <- read_csv("data/gwf_fa_strength.csv") 
d2 <- read_csv("data/gwf_irt.csv") %>% 
  rename(gwf_casename = case) %>%
  dplyr::select(gwf_casename, transition, gwf_irt) 

gwf <- read_csv("data/panel.csv") %>% 
  left_join(d1, by = c("gwf_casename", "year")) %>% 
  left_join(d2, by = c("gwf_casename")) %>% 
  group_by(gwf_casename) %>% 
  mutate_at(.vars = vars(gwf_incumbent_strength), ~na.locf0(.)) %>% 
  group_by(countryname) %>% 
  mutate_at(.vars = vars(lag_milexp, lgupop, e_migdppcln, e_regionpol, e_total_resources_income_pc, islam50, e_peaveduc, ginibestestimateextrapolatednly, LagRuralInequality, larea, gwf_priorduration, e_coups, EFindex, theta_mean), ~na.approx(., maxgap = 5, rule = 2)) %>% 
  mutate_at(.vars = vars(lag_milexp, lgupop, e_migdppcln, e_regionpol, e_total_resources_income_pc, islam50, e_peaveduc, ginibestestimateextrapolatednly, LagRuralInequality, larea, gwf_priorduration, e_coups, EFindex, theta_mean), ~na.locf0(.)) 

gwf$v2pssunpar <- scales::rescale(gwf$v2pssunpar, to = c(0, 1), from = range(gwf$v2pssunpar, na.rm = TRUE, finite = TRUE))

gwf %<>%
  mutate(v2pssunpar = ifelse(is.na(v2xps_party), 0, v2pssunpar)) %>% 
  mutate_if(is.numeric, ~replace(., is.nan(.), NA)) %>% 
  mutate(coups = cumsum(e_coups)) %>% 
  dplyr::select(countryname, year, gwf_casename, gwf_regimetype, v2xps_party, v2pssunpar, gwf_incumbent_strength, 
                gwf_irt, transition,
                e_regionpol, e_migdppcln, lgmilex_cap, lag_milexp, e_peaveduc, lgupop,
                ginibestestimateextrapolatednly, LagRuralInequality, e_total_resources_income_pc, islam50, larea, gwf_priorduration, coups, gwf_number, EFindex, theta_mean) %>%
  mutate(gwf_incumbent_strength = ifelse(is.na(gwf_incumbent_strength), v2pssunpar, gwf_incumbent_strength)) %>% 
  rename(incumbent_strength = gwf_incumbent_strength) %>% 
  filter(transition == year, gwf_casename %notin% c("Bolivia 43-46", "Thailand 57-73", "Pakistan 58-71", "Haiti 88-90",
                                                    "Sudan 58-64", "Ecuador 63-66", "Venezuela 48-58")) %>% 
  ungroup() %>% 
  unique() 

lm(gwf_irt ~ incumbent_strength + ginibestestimateextrapolatednly + LagRuralInequality + lag_milexp + larea + lgupop + e_migdppcln + e_total_resources_income_pc + islam50 + e_regionpol + e_peaveduc, data = gwf) %>% 
  prediction(.) %>% 
  select(incumbent_strength, fitted, se.fitted) %>% 
  ggplot(., aes(x = incumbent_strength, y = fitted)) +
  geom_smooth(color = "black") +
  geom_rug(sides = "b") +
  theme_bw() +
  ylab("Bounded Democracy") +
  ylim(0,1) +
  xlim(0,1) +
  xlab("Party Strength") +
  theme(legend.position="bottom") 

ggsave(filename = "figures/appendix_a_fig1_2.jpg", dpi = 500, width = 5, height = 5) 


#PAR predicted plots Figure 2 & 3#####
#institutionalization#
d1 <- read_csv("data/svolik_fa_pi.csv") 
d2 <- read_csv("data/svolik_irt.csv") %>% 
  rename(svolik_case = case) %>%
  select(svolik_case, transition, svolik_irt)

svolik <- read_csv("data/panel.csv") %>% 
  left_join(d1, by = c("svolik_case", "year")) %>% 
  left_join(d2, by = c("svolik_case")) %>% 
  group_by(svolik_case) %>% 
  mutate_at(.vars = vars(svolik_incumbent_pi), ~na.locf0(.)) %>% 
  group_by(countryname) %>% 
  mutate_at(.vars = vars(lag_milexp, lgupop, e_migdppcln, e_regionpol, e_total_resources_income_pc, islam50, e_peaveduc, ginibestestimateextrapolatednly, LagRuralInequality, larea, gwf_priorduration, e_coups, EFindex, theta_mean), ~na.approx(., maxgap = 5, rule = 2)) %>% 
  mutate_at(.vars = vars(lag_milexp, lgupop, e_migdppcln, e_regionpol, e_total_resources_income_pc, islam50, e_peaveduc, ginibestestimateextrapolatednly, LagRuralInequality, larea, gwf_priorduration, e_coups, EFindex, theta_mean), ~na.locf0(.)) %>% 
  mutate(v2xps_party = ifelse(is.na(v2xps_party), 0, v2xps_party)) %>% 
  mutate_if(is.numeric, ~replace(., is.nan(.), NA)) %>% 
  mutate(coups = cumsum(e_coups)) %>% 
  select(countryname, year, svolik_case, svolik_regime, v2xps_party, v2pssunpar, svolik_incumbent_pi, svolik_irt, transition,
         e_regionpol, e_migdppcln, lgmilex_cap, lag_milexp, e_peaveduc, lgupop,
         ginibestestimateextrapolatednly, LagRuralInequality, e_total_resources_income_pc, islam50, larea,
         svolik_priorduration, coups, svolik_number, EFindex, theta_mean) %>%
  mutate(svolik_incumbent_pi = ifelse(is.na(svolik_incumbent_pi), v2xps_party, svolik_incumbent_pi)) %>% 
  rename(incumbent_pi = svolik_incumbent_pi)  %>% 
  filter(transition == year, svolik_case %notin% c("Bolivia 1946-1946", "Thailand 1971-1972", "Pakistan 1958-1972", "Haiti 1986-1990", 
                                                   "Sudan 1959-1964", "Ecuador 1963-1965", "Haiti 1950-1956",
                                                   "Comoros 1999-2001", "Venezuela 1948-1949", "Burkina Faso (Upper Volta) 1966-1971", "Peru 1949-1956"))%>% 
  ungroup() %>% 
  unique()

lm(svolik_irt ~ incumbent_pi + ginibestestimateextrapolatednly + LagRuralInequality + lag_milexp + larea + lgupop + e_migdppcln + e_total_resources_income_pc + islam50 + e_regionpol + e_peaveduc + svolik_priorduration + coups + svolik_number + EFindex + theta_mean, data = svolik) %>% 
  prediction(.) %>% 
  select(incumbent_pi, fitted, se.fitted) %>% 
  ggplot(., aes(x = incumbent_pi, y = fitted)) +
  geom_smooth(color = "black") +
  geom_rug(sides = "b") +
  theme_bw() +
  ylab("Bounded Democracy") +
  ylim(0,1) +
  xlim(0,1) +
  xlab("Party Institutionalization") +
  theme(legend.position="bottom") 

ggsave(filename = "figures/appendix_a_fig2_1.jpg", dpi = 500, width = 5, height = 5) 

#strength models ##
d1 <- read_csv("data/svolik_fa_strength.csv") 
d2 <- read_csv("data/svolik_irt.csv") %>% 
  rename(svolik_case = case) %>%
  dplyr::select(svolik_case, transition, svolik_irt)

svolik <- read_csv("data/panel.csv") %>% 
  left_join(d1, by = c("svolik_case", "year")) %>% 
  left_join(d2, by = c("svolik_case")) %>% 
  group_by(svolik_case) %>% 
  mutate_at(.vars = vars(svolik_incumbent_strength), ~na.locf0(.)) %>% 
  group_by(countryname) %>% 
  mutate_at(.vars = vars(lag_milexp, lgupop, e_migdppcln, e_regionpol, e_total_resources_income_pc, islam50, e_peaveduc, ginibestestimateextrapolatednly, LagRuralInequality, larea, gwf_priorduration, e_coups, EFindex, theta_mean), ~na.approx(., maxgap = 5, rule = 2)) %>% 
  mutate_at(.vars = vars(lag_milexp, lgupop, e_migdppcln, e_regionpol, e_total_resources_income_pc, islam50, e_peaveduc, ginibestestimateextrapolatednly, LagRuralInequality, larea, gwf_priorduration, e_coups, EFindex, theta_mean), ~na.locf0(.)) 

svolik$v2pssunpar <- scales::rescale(svolik$v2pssunpar, to = c(0, 1), from = range(svolik$v2pssunpar, na.rm = TRUE, finite = TRUE))

svolik %<>%  
  mutate(v2pssunpar = ifelse(is.na(v2xps_party), 0, v2pssunpar))%>% 
  mutate_if(is.numeric, ~replace(., is.nan(.), NA)) %>% 
  mutate(coups = cumsum(e_coups)) %>% 
  dplyr::select(countryname, year, svolik_case, svolik_regime, v2xps_party, v2pssunpar, svolik_incumbent_strength, 
                svolik_irt, transition,
                e_regionpol, e_migdppcln, lgmilex_cap, lag_milexp, e_peaveduc, lgupop,
                ginibestestimateextrapolatednly, LagRuralInequality, e_total_resources_income_pc, islam50, larea,
                svolik_priorduration, coups, svolik_number, EFindex, theta_mean) %>%
  mutate(svolik_incumbent_strength = ifelse(is.na(svolik_incumbent_strength), v2pssunpar, svolik_incumbent_strength)) %>% 
  rename(incumbent_strength = svolik_incumbent_strength) %>% 
  filter(transition == year, svolik_case %notin% c("Bolivia 1946-1946", "Thailand 1971-1972", "Pakistan 1958-1972", "Haiti 1986-1990", 
                                                   "Sudan 1959-1964", "Ecuador 1963-1965", "Haiti 1950-1956",
                                                   "Comoros 1999-2001", "Venezuela 1948-1949", "Burkina Faso (Upper Volta) 1966-1971", "Peru 1949-1956"))%>% 
  ungroup() %>% 
  unique()

lm(svolik_irt ~ incumbent_strength + ginibestestimateextrapolatednly + LagRuralInequality + lag_milexp + larea + lgupop + e_migdppcln + e_total_resources_income_pc + islam50 + e_regionpol + e_peaveduc + svolik_priorduration + coups + svolik_number + EFindex + theta_mean, data = svolik) %>% 
  prediction(.) %>% 
  select(incumbent_strength, fitted, se.fitted, svolik_case) %>% 
  ggplot(., aes(x = incumbent_strength, y = fitted, label = svolik_case)) +
  #geom_label() +
  geom_smooth(color = "black") +
  geom_rug(sides = "b") +
  theme_bw() +
  ylab("Bounded Democracy") +
  ylim(0,1) +
  xlim(0,1) +
  xlab("Party Strength") +
  theme(legend.position="bottom") 

ggsave(filename = "figures/appendix_a_fig2_2.jpg", dpi = 500, width = 5, height = 5) 

#remove Central African Republic Figure 3
svolik_nocar <- filter(svolik, svolik_case != "Central African Republic 1986-1993")

lm(svolik_irt ~ incumbent_strength + ginibestestimateextrapolatednly + LagRuralInequality + lag_milexp + larea + lgupop + e_migdppcln + e_total_resources_income_pc + islam50 + e_regionpol + e_peaveduc + svolik_priorduration + coups + svolik_number + EFindex + theta_mean, data = svolik_nocar) %>% 
  prediction(.) %>% 
  select(incumbent_strength, fitted, se.fitted, svolik_case) %>% 
  ggplot(., aes(x = incumbent_strength, y = fitted, label = svolik_case)) +
  #geom_label() +
  geom_smooth(color = "black") +
  geom_rug(sides = "b") +
  theme_bw() +
  ylab("Bounded Democracy") +
  ylim(0,1) +
  xlim(0,1) +
  xlab("Party Strength") +
  theme(legend.position="bottom") 

ggsave(filename = "figures/appendix_a_fig3.jpg", dpi = 500, width = 5, height = 5) 


#WTH predicted plots Figure 4####
d1 <- read_csv("data/wth_fa_pi.csv") 
d2 <- read_csv("data/wth_irt.csv") %>% 
  rename(wth_case = case) %>%
  select(wth_case, transition, wth_irt)

wth <- read_csv("data/panel.csv") %>% 
  left_join(d1, by = c("wth_case", "year")) %>% 
  left_join(d2, by = c("wth_case")) %>% 
  group_by(wth_case) %>% 
  mutate_at(.vars = vars(wth_incumbent_pi), ~na.locf0(.)) %>% 
  group_by(countryname) %>% 
  mutate_at(.vars = vars(lag_milexp, lgupop, e_migdppcln, e_regionpol, e_total_resources_income_pc, islam50, e_peaveduc, ginibestestimateextrapolatednly, LagRuralInequality, larea, gwf_priorduration, e_coups, EFindex, theta_mean), ~na.approx(., maxgap = 5, rule = 2)) %>% 
  mutate_at(.vars = vars(lag_milexp, lgupop, e_migdppcln, e_regionpol, e_total_resources_income_pc, islam50, e_peaveduc, ginibestestimateextrapolatednly, LagRuralInequality, larea, gwf_priorduration, e_coups, EFindex, theta_mean), ~na.locf0(.)) %>% 
  mutate(v2xps_party = ifelse(is.na(v2xps_party), 0, v2xps_party)) %>% 
  mutate_if(is.numeric, ~replace(., is.nan(.), NA)) %>% 
  mutate(coups = cumsum(e_coups)) %>% 
  select(countryname, year, wth_case, wth_regime, v2xps_party, v2pssunpar, wth_incumbent_pi, wth_irt, transition,
         e_regionpol, e_migdppcln, lgmilex_cap, lag_milexp, e_peaveduc, lgupop,
         ginibestestimateextrapolatednly, LagRuralInequality, e_total_resources_income_pc, islam50, larea, 
         wth_priorduration, coups, wth_number, EFindex, theta_mean) %>%
  mutate(wth_incumbent_pi = ifelse(is.na(wth_incumbent_pi), v2xps_party, wth_incumbent_pi)) %>% 
  rename(incumbent_pi = wth_incumbent_pi)  %>% 
  filter(transition == year, wth_case %notin% c("Haiti 1986-1989", "Comoros 1999-2001")) %>% 
  ungroup() %>% 
  unique()

lm(wth_irt ~ incumbent_pi + ginibestestimateextrapolatednly + LagRuralInequality + lag_milexp + larea + lgupop + e_migdppcln + e_total_resources_income_pc + islam50 + e_regionpol + e_peaveduc + wth_priorduration + coups + wth_number + EFindex + theta_mean, data = wth) %>% 
  prediction(.) %>% 
  select(incumbent_pi, fitted, se.fitted) %>% 
  ggplot(., aes(x = incumbent_pi, y = fitted)) +
  geom_smooth(color = "black") +
  geom_rug(sides = "b") +
  theme_bw() +
  ylab("Bounded Democracy") +
  ylim(0,1) +
  xlim(0,1) +
  xlab("Party Institutionalization") +
  theme(legend.position="bottom") 

ggsave(filename = "figures/appendix_a_fig4_1.jpg", dpi = 500, width = 5, height = 5) 

#strength models #
d1 <- read_csv("data/wth_fa_strength.csv") 
d2 <- read_csv("data/wth_irt.csv") %>% 
  rename(wth_case = case) %>%
  dplyr::select(wth_case, transition, wth_irt)

wth <- read_csv("data/panel.csv") %>% 
  left_join(d1, by = c("wth_case", "year")) %>% 
  left_join(d2, by = c("wth_case")) %>% 
  group_by(wth_case) %>% 
  mutate_at(.vars = vars(wth_incumbent_strength), ~na.locf0(.)) %>% 
  group_by(countryname) %>% 
  mutate_at(.vars = vars(lag_milexp, lgupop, e_migdppcln, e_regionpol, e_total_resources_income_pc, islam50, e_peaveduc, ginibestestimateextrapolatednly, LagRuralInequality, larea, gwf_priorduration, e_coups, EFindex, theta_mean), ~na.approx(., maxgap = 5, rule = 2)) %>% 
  mutate_at(.vars = vars(lag_milexp, lgupop, e_migdppcln, e_regionpol, e_total_resources_income_pc, islam50, e_peaveduc, ginibestestimateextrapolatednly, LagRuralInequality, larea, gwf_priorduration, e_coups, EFindex, theta_mean), ~na.locf0(.)) 

wth$v2pssunpar <- scales::rescale(wth$v2pssunpar, to = c(0, 1), from = range(wth$v2pssunpar, na.rm = TRUE, finite = TRUE))

wth %<>%
  mutate(v2pssunpar = ifelse(is.na(v2xps_party), 0, v2pssunpar)) %>% 
  mutate_if(is.numeric, ~replace(., is.nan(.), NA)) %>% 
  mutate(coups = cumsum(e_coups)) %>% 
  dplyr::select(countryname, year, wth_case, wth_regime, v2xps_party, v2pssunpar, wth_incumbent_strength, 
                wth_irt, transition,
                e_regionpol, e_migdppcln, lgmilex_cap, lag_milexp, e_peaveduc, lgupop,
                ginibestestimateextrapolatednly, LagRuralInequality, e_total_resources_income_pc, islam50, larea, 
                wth_priorduration, coups, wth_number, EFindex, theta_mean) %>%
  mutate(wth_incumbent_strength = ifelse(is.na(wth_incumbent_strength), v2pssunpar, wth_incumbent_strength)) %>% 
  rename(incumbent_strength = wth_incumbent_strength) %>% 
  filter(transition == year, wth_case %notin% c("Haiti 1986-1989", "Comoros 1999-2001")) %>% 
  ungroup() %>% 
  unique() 

lm(wth_irt ~ incumbent_strength + ginibestestimateextrapolatednly + LagRuralInequality + lag_milexp + larea + lgupop + e_migdppcln + e_total_resources_income_pc + islam50 + e_regionpol + e_peaveduc + wth_priorduration + coups + wth_number + EFindex + theta_mean, data = wth) %>% 
  prediction(.) %>% 
  select(incumbent_strength, fitted, se.fitted, wth_case) %>% 
  ggplot(., aes(x = incumbent_strength, y = fitted, label = wth_case)) +
  geom_smooth(color = "black") +
  geom_rug(sides = "b") +
  theme_bw() +
  ylab("Bounded Democracy") +
  ylim(0,1) +
  xlim(0,1) +
  xlab("Party Strength") +
  theme(legend.position="bottom") 

ggsave(filename = "figures/appendix_a_fig4_2.jpg", dpi = 500, width = 5, height = 5) 

#PAR corporate specific models - Figures 5 & 6####
#party institutionalization 
d1 <- read_csv("data/svolik_fa_pi.csv") 
d2 <- read_csv("data/svolik_irt.csv") %>% 
  rename(svolik_case = case) %>%
  select(svolik_case, transition, svolik_irt)

svolik <- read_csv("data/panel.csv") %>% 
  left_join(d1, by = c("svolik_case", "year")) %>% 
  left_join(d2, by = c("svolik_case")) %>% 
  group_by(svolik_case) %>% 
  mutate_at(.vars = vars(svolik_incumbent_pi), ~na.locf0(.)) %>% 
  group_by(countryname) %>% 
  mutate_at(.vars = vars(lag_milexp, lgupop, e_migdppcln, e_regionpol, e_total_resources_income_pc, islam50, e_peaveduc, ginibestestimateextrapolatednly, LagRuralInequality, larea, gwf_priorduration, e_coups, EFindex, theta_mean), ~na.approx(., maxgap = 5, rule = 2)) %>% 
  mutate_at(.vars = vars(lag_milexp, lgupop, e_migdppcln, e_regionpol, e_total_resources_income_pc, islam50, e_peaveduc, ginibestestimateextrapolatednly, LagRuralInequality, larea, gwf_priorduration, e_coups, EFindex, theta_mean), ~na.locf0(.)) %>% 
  mutate(v2xps_party = ifelse(is.na(v2xps_party), 0, v2xps_party)) %>% 
  mutate_if(is.numeric, ~replace(., is.nan(.), NA)) %>% 
  mutate(coups = cumsum(e_coups)) %>% 
  select(countryname, year, svolik_case, svolik_regime, v2xps_party, v2pssunpar, svolik_incumbent_pi, svolik_irt, transition,
         e_regionpol, e_migdppcln, lgmilex_cap, lag_milexp, e_peaveduc, lgupop,
         ginibestestimateextrapolatednly, LagRuralInequality, e_total_resources_income_pc, islam50, larea,
         svolik_priorduration, coups, svolik_number, EFindex, theta_mean) %>%
  mutate(svolik_incumbent_pi = ifelse(is.na(svolik_incumbent_pi), v2xps_party, svolik_incumbent_pi)) %>% 
  rename(incumbent_pi = svolik_incumbent_pi)  %>% 
  filter(transition == year, svolik_case %notin% c("Bolivia 1946-1946", "Thailand 1971-1972", "Pakistan 1958-1972", "Haiti 1986-1990", 
                                                   "Sudan 1959-1964", "Ecuador 1963-1965", "Haiti 1950-1956",
                                                   "Comoros 1999-2001", "Venezuela 1948-1949", "Burkina Faso (Upper Volta) 1966-1971", "Peru 1949-1956"))%>% 
  ungroup() %>% 
  unique()

d3 <- read_csv("data/Svolik clean.csv") %>% 
  dplyr::select(countryname, year, par_military, par_regime) %>% 
  mutate(par_military = ifelse(countryname == "Brazil" & year %in% 1965:1986, "corporate", par_military)) %>% 
  group_by(countryname) %>% 
  mutate(par_military = ifelse(is.na(par_military), lag(par_military), par_military)) %>% 
  mutate(svolik_regime = ifelse(par_regime == "democracy", "democracy", par_military),
         svolik_regime = ifelse(par_regime == "no authority" & is.na(svolik_regime), par_regime, svolik_regime), 
         svolik_regime = ifelse(is.na(svolik_regime), par_regime, svolik_regime))  %>%   
  mutate(svolik_number = ifelse(lag(svolik_regime) != svolik_regime, 1, 0), svolik_number = ifelse(is.na(svolik_number), 0, svolik_number),
         svolik_number = cumsum(svolik_number)+1) %>% 
  mutate(svolik_regime = ifelse(svolik_regime %in% c("corporate", "personal", "indirect"), "military", svolik_regime), 
         svolik_regime = ifelse(is.na(svolik_regime) & par_regime == "dictatorship", "dictatorship", svolik_regime)) %>% 
  dplyr::select(countryname, year, svolik_regime, svolik_number, par_military) %>% 
  group_by(countryname, svolik_number) %>% 
  mutate(eyear = ifelse(svolik_number != lead(svolik_number), year, NA),
         syear = ifelse(svolik_number != lag(svolik_number), year, NA),
         syear = na.locf0(syear), syear = ifelse(is.na(syear), first(year), syear),
         eyear = na.locf0(eyear), eyear = ifelse(is.na(eyear), last(year), eyear)) %>% 
  group_by(countryname, svolik_number) %>% 
  mutate(svolik_case = paste(syear, eyear, sep = "-"),
         svolik_case = paste(countryname, svolik_case))  %>% 
  dplyr::select(-year) %>% 
  unique() %>% 
  group_by(countryname) %>% 
  mutate(svolik_next = lead(svolik_regime)) %>% 
  ungroup() %>% 
  dplyr::select(svolik_case, par_military)


df <- left_join(svolik, d3, by = c("svolik_case")) %>% 
  mutate(par_military = ifelse(svolik_case == "Korea, Republic of 1961-1988", "corporate", par_military)) %>% 
  filter(par_military == "corporate")

m1 <- lm_robust(svolik_irt ~ incumbent_pi, data = df, se_type = "HC1") %>% 
  tidy(.) %>% 
  mutate(Model = "Bivariate") 

m2 <- lm_robust(svolik_irt ~ incumbent_pi + lag_milexp + larea + lgupop + e_migdppcln + e_regionpol + e_total_resources_income_pc + islam50 + e_peaveduc + svolik_priorduration + coups + svolik_number + EFindex + theta_mean, data = df, se_type = "HC1") %>% 
  tidy(.) %>% 
  mutate(Model = "Base") 

m3 <- lm_robust(svolik_irt ~ incumbent_pi + ginibestestimateextrapolatednly + LagRuralInequality + lag_milexp + larea + lgupop + e_migdppcln + e_total_resources_income_pc + islam50 + e_regionpol + e_peaveduc + svolik_priorduration + coups + svolik_number + EFindex + theta_mean, data = df, se_type = "HC1") %>% 
  tidy(.) %>% 
  mutate(Model = "Inequality")

#coefficient plot for institutionalization figure 5
rbind(m1, m2, m3) %>% 
  filter(term != "(Intercept)") %>% 
  mutate(conf.low = round(conf.low, 3), conf.high = round(conf.high, 3)) %>% 
  mutate(term = factor(term, levels = c("theta_mean", "EFindex", "svolik_number", "coups", "svolik_priorduration", "lag_milexp", "larea", "lgupop", "e_migdppcln", "e_regionpol", "e_total_resources_income_pc", "islam50", 
                                        "e_peaveduc", "ginibestestimateextrapolatednly", "LagRuralInequality", "incumbent_pi"))) %>% 
  mutate(Model = factor(Model, levels = c("Bivariate", "Base", "Inequality"))) %>% 
  mutate(term = recode(term, "theta_mean" = "Repression", "lag_milexp" = "Military Spending", "larea" = "Area", "lgupop" = "Urban", "e_migdppcln" = "GDP PC", "e_regionpol" = "Region", 
                       "e_total_resources_income_pc" = "Resources", "islam50" = "Islam", "e_peaveduc" = "Education", "ginibestestimateextrapolatednly" = "GINI", 
                       "LagRuralInequality" = "Rural", "incumbent_pi" = "Institutionalization", "svolik_priorduration" = "Duration", "coups" = "Coups", "svolik_number" = "Breakdowns", "EFindex" = "Ethnic Frac")) %>% 
  ggplot(., aes(x = term, y = estimate, ymin = conf.low, ymax = conf.high, color = factor(Model), shape = Model)) +
  geom_pointrange(show.legend = F, position = position_jitterdodge(dodge.width = 0.75, jitter.width = 0.2), size = 0.5) +
  geom_hline(yintercept=0, linetype = "dashed") +
  geom_line(aes(y = estimate-100)) +
  geom_point(aes(y = estimate-100)) +
  theme_bw() +
  xlab("Party Institutionalization") +
  ylab("Bounded Democracy") +
  ylim(-2, 2) +
  scale_color_grey(name = "Model", start = 0.8, end = 0) +
  theme(legend.background = element_rect(fill = alpha(.1)), legend.position = "bottom") +
  coord_flip() 

ggsave(filename = "figures/appendix_a_fig5_1.jpg", dpi = 500, width = 5, height = 5) 

#predicted plot - institutionalization Figure 6
lm(svolik_irt ~ incumbent_pi + ginibestestimateextrapolatednly + LagRuralInequality + lag_milexp + larea + lgupop + e_migdppcln + e_total_resources_income_pc + islam50 + e_regionpol + e_peaveduc + svolik_priorduration + coups + svolik_number + EFindex + theta_mean, data = df) %>% 
  prediction(.) %>% 
  select(incumbent_pi, fitted, se.fitted, svolik_case) %>% 
  ggplot(., aes(x = incumbent_pi, y = fitted, label = svolik_case)) +
  geom_smooth(color = "black") +
  geom_rug(sides = "b") +
  theme_bw() +
  ylab("Bounded Democracy") +
  ylim(0,1) +
  xlim(0,1) +
  xlab("Party Institutionalization") +
  theme(legend.position="bottom") 

ggsave(filename = "figures/appendix_a_fig6_1.jpg", dpi = 500, width = 5, height = 5) 

#party strength
d1 <- read_csv("data/svolik_fa_strength.csv") 
d2 <- read_csv("data/svolik_irt.csv") %>% 
  rename(svolik_case = case) %>%
  dplyr::select(svolik_case, transition, svolik_irt)

svolik <- read_csv("data/panel.csv") %>% 
  left_join(d1, by = c("svolik_case", "year")) %>% 
  left_join(d2, by = c("svolik_case")) %>% 
  group_by(svolik_case) %>% 
  mutate_at(.vars = vars(svolik_incumbent_strength), ~na.locf0(.)) %>% 
  group_by(countryname) %>% 
  mutate_at(.vars = vars(lag_milexp, lgupop, e_migdppcln, e_regionpol, e_total_resources_income_pc, islam50, e_peaveduc, ginibestestimateextrapolatednly, LagRuralInequality, larea, gwf_priorduration, e_coups, EFindex, theta_mean), ~na.approx(., maxgap = 5, rule = 2)) %>% 
  mutate_at(.vars = vars(lag_milexp, lgupop, e_migdppcln, e_regionpol, e_total_resources_income_pc, islam50, e_peaveduc, ginibestestimateextrapolatednly, LagRuralInequality, larea, gwf_priorduration, e_coups, EFindex, theta_mean), ~na.locf0(.)) 

svolik$v2pssunpar <- scales::rescale(svolik$v2pssunpar, to = c(0, 1), from = range(svolik$v2pssunpar, na.rm = TRUE, finite = TRUE))

svolik %<>%  
  mutate(v2pssunpar = ifelse(is.na(v2xps_party), 0, v2pssunpar))%>% 
  mutate_if(is.numeric, ~replace(., is.nan(.), NA)) %>% 
  mutate(coups = cumsum(e_coups)) %>% 
  dplyr::select(countryname, year, svolik_case, svolik_regime, v2xps_party, v2pssunpar, svolik_incumbent_strength, 
                svolik_irt, transition,
                e_regionpol, e_migdppcln, lgmilex_cap, lag_milexp, e_peaveduc, lgupop,
                ginibestestimateextrapolatednly, LagRuralInequality, e_total_resources_income_pc, islam50, larea,
                svolik_priorduration, coups, svolik_number, EFindex, theta_mean) %>%
  mutate(svolik_incumbent_strength = ifelse(is.na(svolik_incumbent_strength), v2pssunpar, svolik_incumbent_strength)) %>% 
  rename(incumbent_strength = svolik_incumbent_strength) %>% 
  filter(transition == year, svolik_case %notin% c("Bolivia 1946-1946", "Thailand 1971-1972", "Pakistan 1958-1972", "Haiti 1986-1990", 
                                                   "Sudan 1959-1964", "Ecuador 1963-1965", "Haiti 1950-1956",
                                                   "Comoros 1999-2001", "Venezuela 1948-1949", "Burkina Faso (Upper Volta) 1966-1971", "Peru 1949-1956"))%>% 
  ungroup() %>% 
  unique()

df <- left_join(svolik, d3, by = c("svolik_case")) %>% 
  mutate(par_military = ifelse(svolik_case == "Korea, Republic of 1961-1988", "corporate", par_military)) %>% 
  filter(par_military == "corporate")

#coefficient plot - strength Figure 5
m1 <- lm_robust(svolik_irt ~ incumbent_strength, data = df, se_type = "HC1") %>% 
  tidy(.) %>% 
  mutate(Model = "Bivariate") 

m2 <- lm_robust(svolik_irt ~ incumbent_strength + lag_milexp + larea + lgupop + e_migdppcln + e_regionpol + e_total_resources_income_pc + islam50 + e_peaveduc + svolik_priorduration + coups + svolik_number + EFindex + theta_mean, data = df, se_type = "HC1") %>% 
  tidy(.) %>% 
  mutate(Model = "Base") 

m3 <- lm_robust(svolik_irt ~ incumbent_strength + ginibestestimateextrapolatednly + LagRuralInequality + lag_milexp + larea + lgupop + e_migdppcln + e_total_resources_income_pc + islam50 + e_regionpol + e_peaveduc + svolik_priorduration + coups + svolik_number + EFindex + theta_mean, data = df, se_type = "HC1") %>% 
  tidy(.) %>% 
  mutate(Model = "Inequality")

rbind(m1, m2, m3) %>% 
  filter(term != "(Intercept)") %>% 
  mutate(conf.low = round(conf.low, 3), conf.high = round(conf.high, 3)) %>% 
  mutate(term = factor(term, levels = c("theta_mean", "EFindex", "svolik_number", "coups", "svolik_priorduration", "lag_milexp", "larea", "lgupop", "e_migdppcln", "e_regionpol", "e_total_resources_income_pc", "islam50", 
                                        "e_peaveduc", "ginibestestimateextrapolatednly", "LagRuralInequality", "incumbent_strength"))) %>% 
  mutate(Model = factor(Model, levels = c("Bivariate", "Base", "Inequality"))) %>% 
  mutate(term = recode(term, "theta_mean" = "Repression", "lag_milexp" = "Military Spending", "larea" = "Area", "lgupop" = "Urban", "e_migdppcln" = "GDP PC", "e_regionpol" = "Region", 
                       "e_total_resources_income_pc" = "Resources", "islam50" = "Islam", "e_peaveduc" = "Education", "ginibestestimateextrapolatednly" = "GINI", 
                       "LagRuralInequality" = "Rural", "incumbent_strength" = "Strength", "svolik_priorduration" = "Duration", "coups" = "Coups", "svolik_number" = "Breakdowns", "EFindex" = "Ethnic Frac")) %>% 
  ggplot(., aes(x = term, y = estimate, ymin = conf.low, ymax = conf.high, color = factor(Model), shape = Model)) +
  geom_pointrange(show.legend = F, position = position_jitterdodge(dodge.width = 0.75, jitter.width = 0.2), size = 0.5) +
  geom_hline(yintercept=0, linetype = "dashed") +
  geom_line(aes(y = estimate-100)) +
  geom_point(aes(y = estimate-100)) +
  theme_bw() +
  xlab("Party Strength") +
  ylab("Bounded Democracy") +
  ylim(-2, 2) +
  scale_color_grey(name = "Model", start = 0.8, end = 0) +
  theme(legend.background = element_rect(fill = alpha(.1)), legend.position = "bottom") +
  coord_flip() 

ggsave(filename = "figures/appendix_a_fig5_2.jpg", dpi = 500, width = 5, height = 5) 

#predicted plot - strength Figure 6
lm(svolik_irt ~ incumbent_strength + ginibestestimateextrapolatednly + LagRuralInequality + lag_milexp + larea + lgupop + e_migdppcln + e_total_resources_income_pc + islam50 + e_regionpol + e_peaveduc + svolik_priorduration + coups + svolik_number + EFindex + theta_mean, data = df) %>% 
  prediction(.) %>% 
  select(incumbent_strength, fitted, se.fitted, svolik_case) %>% 
  ggplot(., aes(x = incumbent_strength, y = fitted, label = svolik_case)) +
  geom_smooth(color = "black") +
  geom_rug(sides = "b") +
  theme_bw() +
  ylab("Bounded Democracy") +
  ylim(0,1) +
  xlim(0,1) +
  xlab("Party Strength") +
  theme(legend.position="bottom") 

ggsave(filename = "figures/appendix_a_fig6_2.jpg", dpi = 500, width = 5, height = 5) 

#DD plots Figure 7 ####
#institutionalization#
d1 <- read_csv("data/dd_fa_pi.csv") 
d2 <- read_csv("data/dd_irt.csv") %>% 
  rename(dd_case = case) %>%
  select(dd_case, transition, dd_irt)

dd <- read_csv("data/panel.csv") %>% 
  left_join(d1, by = c("dd_case", "year")) %>% 
  left_join(d2, by = c("dd_case")) %>% 
  group_by(dd_case) %>% 
  mutate_at(.vars = vars(dd_incumbent_pi), ~na.locf0(.)) %>% 
  group_by(countryname) %>% 
  mutate_at(.vars = vars(lag_milexp, lgupop, e_migdppcln, e_regionpol, e_total_resources_income_pc, islam50, e_peaveduc, ginibestestimateextrapolatednly, LagRuralInequality, larea, gwf_priorduration, e_coups, EFindex, theta_mean), ~na.approx(., maxgap = 5, rule = 2)) %>% 
  mutate_at(.vars = vars(lag_milexp, lgupop, e_migdppcln, e_regionpol, e_total_resources_income_pc, islam50, e_peaveduc, ginibestestimateextrapolatednly, LagRuralInequality, larea, gwf_priorduration, e_coups, EFindex, theta_mean), ~na.locf0(.)) %>% 
  mutate(v2xps_party = ifelse(is.na(v2xps_party), 0, v2xps_party)) %>% 
  mutate_if(is.numeric, ~replace(., is.nan(.), NA)) %>% 
  mutate(coups = cumsum(e_coups)) %>% 
  select(countryname, year, dd_case, dd_regime, v2xps_party, v2pssunpar, dd_incumbent_pi, dd_irt, transition,
         e_regionpol, e_migdppcln, lgmilex_cap, lag_milexp, e_peaveduc, lgupop,
         ginibestestimateextrapolatednly, LagRuralInequality, e_total_resources_income_pc, islam50, larea, 
         dd_priorduration, coups, dd_number, EFindex, theta_mean) %>%
  mutate(dd_incumbent_pi = ifelse(is.na(dd_incumbent_pi), v2xps_party, dd_incumbent_pi)) %>% 
  rename(incumbent_pi = dd_incumbent_pi)  %>% 
  filter(transition == year, dd_case %notin% c("Bolivia 1946-1946", "Thailand 1946-1972", "Pakistan 1958-1970", "Haiti 1986-1989", 
                                               "Sudan 1958-1963", "Ecuador 1963-1965", "Haiti 1950-1955",
                                               "Comoros 1999-2003", "Venezuela 1948-1958")) %>% 
  ungroup() %>% 
  unique()

lm(dd_irt ~ incumbent_pi + ginibestestimateextrapolatednly + LagRuralInequality + lag_milexp + larea + lgupop + e_migdppcln + e_total_resources_income_pc + islam50 + e_regionpol + e_peaveduc + dd_priorduration + coups + dd_number + EFindex + theta_mean, data = dd) %>% 
  prediction(.) %>% 
  select(incumbent_pi, fitted, se.fitted) %>% 
  ggplot(., aes(x = incumbent_pi, y = fitted)) +
  geom_smooth(color = "black") +
  geom_rug(sides = "b") +
  theme_bw() +
  ylab("Bounded Democracy") +
  ylim(0,1) +
  xlim(0,1) +
  xlab("Party Institutionalization") +
  theme(legend.position="bottom") 

ggsave(filename = "figures/appendix_a_fig7_1.jpg", dpi = 500, width = 5, height = 5) 

#strength models 
d1 <- read_csv("data/dd_fa_strength.csv") 
d2 <- read_csv("data/dd_irt.csv") %>% 
  rename(dd_case = case) %>%
  dplyr::select(dd_case, transition, dd_irt)

dd <- read_csv("data/panel.csv") %>% 
  left_join(d1, by = c("dd_case", "year")) %>% 
  left_join(d2, by = c("dd_case")) %>% 
  group_by(dd_case) %>% 
  mutate_at(.vars = vars(dd_incumbent_strength), ~na.locf0(.)) %>% 
  group_by(countryname) %>% 
  mutate_at(.vars = vars(lag_milexp, lgupop, e_migdppcln, e_regionpol, e_total_resources_income_pc, islam50, e_peaveduc, ginibestestimateextrapolatednly, LagRuralInequality, larea, gwf_priorduration, e_coups, EFindex, theta_mean), ~na.approx(., maxgap = 5, rule = 2)) %>% 
  mutate_at(.vars = vars(lag_milexp, lgupop, e_migdppcln, e_regionpol, e_total_resources_income_pc, islam50, e_peaveduc, ginibestestimateextrapolatednly, LagRuralInequality, larea, gwf_priorduration, e_coups, EFindex, theta_mean), ~na.locf0(.)) 

dd$v2pssunpar <- scales::rescale(dd$v2pssunpar, to = c(0, 1), from = range(dd$v2pssunpar, na.rm = TRUE, finite = TRUE))

dd %<>%
  mutate(v2pssunpar = ifelse(is.na(v2xps_party), 0, v2pssunpar)) %>% 
  mutate_if(is.numeric, ~replace(., is.nan(.), NA)) %>% 
  mutate(coups = cumsum(e_coups)) %>% 
  dplyr::select(countryname, year, dd_case, dd_regime, v2xps_party, v2pssunpar, dd_incumbent_strength, 
                dd_irt, transition,
                e_regionpol, e_migdppcln, lgmilex_cap, lag_milexp, e_peaveduc, lgupop,
                ginibestestimateextrapolatednly, LagRuralInequality, e_total_resources_income_pc, islam50, larea, 
                dd_priorduration, coups, dd_number, EFindex, theta_mean) %>%
  mutate(dd_incumbent_strength = ifelse(is.na(dd_incumbent_strength), v2pssunpar, dd_incumbent_strength)) %>% 
  rename(incumbent_strength = dd_incumbent_strength) %>% 
  filter(transition == year, dd_case %notin% c("Bolivia 1946-1946", "Thailand 1946-1972", "Pakistan 1958-1970", "Haiti 1986-1989", 
                                               "Sudan 1958-1963", "Ecuador 1963-1965", "Haiti 1950-1955",
                                               "Comoros 1999-2003", "Venezuela 1948-1958")) %>% 
  ungroup() %>% 
  unique() 

lm(dd_irt ~ incumbent_strength + ginibestestimateextrapolatednly + LagRuralInequality + lag_milexp + larea + lgupop + e_migdppcln + e_total_resources_income_pc + islam50 + e_regionpol + e_peaveduc + dd_priorduration + coups + dd_number + EFindex + theta_mean, data = dd) %>% 
  prediction(.) %>% 
  select(incumbent_strength, fitted, se.fitted, dd_case) %>% 
  ggplot(., aes(x = incumbent_strength, y = fitted, label = dd_case)) +
  geom_smooth(color = "black") +
  geom_rug(sides = "b") +
  theme_bw() +
  ylab("Bounded Democracy") +
  ylim(0,1) +
  xlim(0,1) +
  xlab("Party Strength") +
  theme(legend.position="bottom") 

ggsave(filename = "figures/appendix_a_fig7_2.jpg", dpi = 500, width = 5, height = 5) 


